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1  Introduction 

The  BUCKDEL  software  package  consists  of  two  programs.  The  first  program 
is  called  BUCKGEN  and  it  is  used  to  generate  the  input  data  file  for  the  main 
program.  The  main  program,  called  BUCKDEL,  performs  geometrically 
nonlinear  analysis  of  stiffened  laminated  composite  plates  with  and  without 
delaminations.  BUCKDEL  allows  the  user  to  perform: 

•  a  linear  static  solution; 

•  a  linear  buckling  (eigenvalue)  analysis;  and 

•  a  nonlinear  post-buckling  analysis  through  both  limit  and  bifurcation  points. 
There  are  two  finite  elements  in  BUCKDEL;  a  quasi-conforming  triangular 
laminated  plate  element  based  on  a  refined  first  order  shear  theory  with  seven 
degrees  of  freedom  per  node  and  a  fiber-reinforced  beam  element  with  seven 
degrees  of  freedom  per  node.  The  behavior  of  all  materials  is  assumed  to  be 
linearly  elastic. 


The  multi-domain  modeling  method  is  used  for  the  analysis  of  laminated  plates 
containing  delaminations.  In  this  model,  the  delaminate,  base  and  the 
undelaminate  are  modeled  as  3  distinct  plates  (see  Figure  1-1).  On  the 
delamination  front,  Mindlin’s  deformation  assumption  is  applied  to  obtain  the 
relationship  between  the  displacements  in  the  delaminate,  base  and 
undelaminate. 


Figure  1-1:  Multi-domain  model  of  a  laminated  plate  containing  a  delamination 


2  Disk  Contents  and  Installation  Procedure 

This  section  describes  the  steps  that  must  be  followed  to  use  the  BUCKDEL 
software  package.  The  disk  accompanying  this  manual  contains: 

1.  executable  versions  of  BUCKGEN  (buckgen.exe)  and  BUCKDEL 

(buckdel.exe); 
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2.  sample  input  files;  and 

3.  the  file  readme.txt. 

The  sample  input  files  are  described  in  Section  7  of  this  manual.  They  can  be 
used  to  verify  that  BUCKDEL  is  functioning  correctly  on  your  system  once  it  is 
installed.  The  file  readme.txt  contains  installation  and  program  execution 
instructions  specific  to  your  machine.  All  the  files  on  the  disk  should  be  loaded 
onto  your  computer.  It  will  be  convenient  to  create  a  separate  directory  for  these 
files.  The  procedure  required  to  load  these  files  onto  your  computer  will  be 
system  dependent. 


3  Running  BUCKDEL 

To  execute  BUCKDEL,  the  user  must  first  create  an  input  data  file.  Therefore, 
the  interactive  input  file  generator,  BUCKGEN,  must  be  run  prior  to  executing 
BUCKDEL.  The  procedures  for  running  BUCKGEN  are  described  in  detail  in 
Section  4.  While  BUCKGEN  automates  the  creation  of  the  BUCKDEL  input 
file,  there  are  several  input  file  parameters  which  must  be  manually  added 
by  the  user.  These  are  discussed  in  Section  4.3.  The  structure  of  the 
BUCKDEL  input  file  created  by  BUCKGEN  is  described  in  detail  in  Section  5. 
Once  an  input  data  file  has  been  prepared,  BUCKDEL  is  easy  to  run.  When 
BUCKDEL  is  executed,  the  user  is  prompted  for  the  name  of  the  input  file. 
BUCKDEL  reads  the  specified  input  file  and  execution  proceeds.  During 
execution,  various  output  files  are  created  by  BUCKDEL.  These  depend  to  some 
extent  on  the  analysis  option  chosen  and  are  described  in  Section  6. 


4  Running  BUCKGEN 

4.1  Introduction 

BUCKGEN  is  an  interactive  pre-processor  that  creates  data  files  that  are  used  as 
input  to  BUCKDEL.  The  following  are  several  general  comments  on  the 
operation  of  BUCKGEN.  It  is  important  that  the  user  be  familiar  with  these 
prior  to  running  BUCKGEN. 

•  To  initiate  BUCKGEN,  the  user  simply  types  buckgen.exe. 

•  The  input  for  BUCKGEN  is  primarily  carried  out  by  stepping  through  a 
series  of  input  menus,  each  containing  several  options.  Once  a  particular 
option  is  selected,  the  program  will  issue  a  number  of  prompts  to  the  user 
requesting  that  data  be  entered  or  that  additional  options  be  selected.  Data 
can  be  easily  modified  or  changed  by  re-selecting  a  particular  option  at  a  later 
stage. 

•  Items  from  a  menu  are  selected  by  entering  the  number  corresponding  to  a 
particular  item. 

•  BUCKGEN  contains  many  features  that  will  assist  the  user  in  generating  an 
input  file  for  BUCKDEL.  For  example,  there  are  a  number  of  warning 
messages  advising  the  user  when  the  input  data  is  inconsistent  with  that 
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expected.  Because  of  its  interactive  nature,  the  program  allows  the  user  to 
re-enter  the  data  when  such  a  warning  message  is  received. 


4.2  BUCKGEN  Main  Menu 

The  main  menu,  which  appears  when  BUCKGEN  is  initiated,  contains  three 
options.  The  user  selects  the  desired  option  by  entering  the  appropriate  number. 
These  options  are  as  follows. 


1  :Control  Data 

2  :Mesh  Generation  Data 

3  :Exit  BUCKGEN 


When  creating  a  data  file,  the  user  must  select  some  or  all  of  these  options. 
Option  1  must  always  be  selected  prior  to  all  other  options  and  this  option 
must  be  completed  before  proceeding  to  the  other  options.  The  remaining 
options  can  be  carried  out  in  a  random  order. 

4.2.1  BUCKGEN  Input  Menu  1:  Control  Data 

Various  control  parameters  required  by  BUCKGEN  and  BUCKDEL  are  entered 
from  Input  Menu  1,  which  appears  when  item  1  is  selected  from  the  Main  Menu. 
Input  Menu  1  contains  four  options.  These  options  are  as  follows. 


1  :Enter  Problem  Title 

2  :Enter  Output  File  Name  (REQUIRED) 

3  :  Enter  Analysis  Option 

4  :  Return  to  Main  Menu 


The  user  must  complete  this  menu  before  proceeding  with  any  other  portions  of 
the  data  file  generation.  Note  that  Option  2  must  be  completed. 

4.2. 1 . 1  Option  1 :  Enter  Problem  Title 

When  this  option  is  selected,  the  user  is  prompted  for  a  problem  title.  This 
option  is  for  the  benefit  of  the  user  and  may  be  skipped.  The  title  must  be  less 
than  sixty  characters. 

4.2. 1 .2  Option  2:  Enter  Output  File  Name 
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The  user  must  always  select  a  name  for  the  output  file.  This  can  have  the  same 
name  as  an  existing  file  or  it  can  be  a  new  file.  It  must  be  less  than  twenty 
characters  in  length. 

4.2. 1.3  Option  3:  Enter  Analysis  Option 

When  this  option  is  selected,  a  new  menu  appears  containing  the  analysis 
choices.  The  analysis  choices  are  as  follows. 


1  :Linear  analysis 

2  :Linear  and  buckling  analysis 

3  :  Buckling  and  post-buckling  analysis 

4  :Nonlinear  analysis  through  both  limit  and  bifurcation  points 

5  :Nonlinear  analysis  through  limit  points  only 


When  a  selection  is  made,  control  automatically  reverts  back  to  Input  Menu  1. 

4.2. 1 .4  Option  4:  Return  to  Main  Menu 

Select  this  option  in  order  to  return  to  the  Main  Menu. 

4.2.2  BUCKGEN  Input  Menu  2:  Mesh  Generation  Data 

The  data  defining  the  finite  element  model  is  entered  through  Input  Menu  2.  The 
data  (i.e.  nodes,  element  connectivities,  material  and  geometric  properties  and 
boundary  conditions)  is  imported  from  a  file  which  is  in  a  NASTRAN  bulk  data 
file  format.  The  options  in  Input  Menu  2  are  as  follows. 


1  :Translate  Finite  Element  Model  Data  File 

2  :View  Nodal  Coordinates 

3  :View  Element  Connectivities 

4  :  Optimize  Node  Numbering 

5  :  Return  to  Main  Menu 


4.2.2. 1  Option  1:  Translate  Finite  Element  Model  Data  File 
BUCKGEN  has  been  set  up  to  read  finite  element  model  data  from  a  NASTRAN 
bulk  data  file  and  translate  it  into  a  form  which  BUCKDEL  can  use.  Thus,  any 
finite  element  pre-processor  which  can  write  a  file  in  this  format  can  be  used  in 
conjunction  with  the  BUCKDEL  package.  The  bulk  data  file  should  contain 
only  information  about  nodes,  element  connectivities,  material  and  geometric 
properties  and  boundary  conditions.  The  file  should  begin  with  “BEGIN  BULK” 
and  end  with  “ENDDATA”.  The  bulk  data  mnemonics  which  are  supported  are: 
GRID,  CTRIA3,  CBAR,  MAT1,  MAT8,  PSHELL,  PBAR,  PCOMP,  FORCE, 


4 


SPC,  and  SPC1.  In  addition,  the  mnemonics  $,  PARAM,  and  CORD  are  read  by 
BUCKGEN  but  data  associated  with  them  is  not  used.  For  proper  translation  to 
take  place,  the  following  guidelines  should  be  observed  when  creating  the  finite 
element  model. 

•  The  finite  element  model  should  be  created  using  a  single  Cartesian 
coordinate  system. 

•  The  finite  element  model  should  be  created  in  the  x-y  plane. 

•  The  plate  should  be  modeled  with  3  noded  triangular  elements. 

•  The  stiffeners  should  be  modeled  with  2  noded  beam  elements. 

•  The  element  connectivities  of  beam  elements  should  be  specified  using  nodes 
which  are  also  used  to  specify  the  connectivities  of  triangular  elements  (i.e. 
there  should  not  be  any  nodes  which  are  referenced  by  beam  elements  only). 

•  The  finite  element  model  should  have  consecutively  numbered  nodes  and 
elements,  each  beginning  with  number  one. 

•  Non-zero  displacement  constraints  are  not  permitted. 

•  Only  nodal  loads  may  be  applied  to  the  model. 

•  An  isotropic  or  orthotropic  material  may  be  specified  for  the  plate. 

•  Ply  angles  are  given  with  respect  to  the  global  x-axis. 

•  Ply  angles  and  thicknesses  are  assumed  to  specified  from  the  “bottom”  or 
z”  face  of  the  plate. 

•  An  isotropic  material  must  be  used  for  the  stiffeners. 

•  Delaminations  are  modeled  by  three  distinct  plates  (see  Figure  1).  For  the 
purposes  of  model  creation,  it  is  most  convenient  to  create  the  undelaminate 
in  the  z=0  plane,  the  base  in  the  z=-a  plane  and  the  delaminate  in  the  z=a 
plane  (a  is  an  arbitrary  number).  Upon  translation,  BUCKGEN  will  map  the 
three  plates  to  the  z=0  plane. 

•  The  undelaminate  should  be  assigned  property  number  1,  the  base  property 
number  2,  and  the  delaminate  property  number  3. 

•  The  number  of  plys  in  the  delaminate  and  base  should  equal  that  of  the 
undelaminate. 

•  The  nodes  of  the  undelaminate,  base  and  delaminate  on  the  delamination 
front  should  have  the  same  x  and  y  coordinates.  After  translation  by 
BUCKGEN  and  mapping  of  the  three  plates  to  the  z-0  plane,  the  nodes  of 
the  three  plates  on  the  delamination  front  will  be  coincident. 

•  The  plate  elements  should  be  symmetric  across  the  delamination  front  if 
energy  release  rate  calculations  are  being  performed  (see  Cards  for  Plate 
Elements  in  Section  5). 

•  Stiffeners  can’t  be  modeled  on  base  and  delaminate  plates.  They  can  only  be 
modeled  on  the  undelaminate  plate. 

When  this  option  is  selected,  a  menu  appears  which  gives  the  user  one  choice 
(translation  from  other  formats  will  be  added  in  the  future). 


1 1  '.Translate  From  NASTRAN  Bulk  Data  File 
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After  making  a  choice,  the  user  is  then  prompted  for  a  file  name.  The  file  name 
must  be  less  than  twenty  characters.  BUCKGEN  will  translate  the  designated 
file  and  return  to  Input  Menu  2.  During  the  translation,  if  data  cards  are  read 
which  BUCKGEN  does  not  recognize,  the  user  will  be  warned. 

4.2.2. 2  Option  2:  View  Nodal  Coordinates 

The  user  can  view  the  nodal  coordinates  of  selected  nodes  by  choosing  this 
option.  Since  the  number  of  nodes  in  the  system  is  often  relatively  large,  the 
user  may  not  wish  to  look  at  all  the  data.  Instead  the  user  can  view  a  range  of 
nodes.  Thus,  when  the  option  is  selected  the  user  is  prompted  for  the  first  and 
last  nodes  in  the  range.  The  coordinates  of  these  nodes  are  then  written  to  the 
screen. 

4.2.2.3  Option  3:  View  Element  Connectivities 

The  user  can  view  the  connectivities  of  selected  elements  by  choosing  this 
option.  Since  the  number  of  elements  in  the  system  is  often  relatively  large,  the 
user  may  not  wish  to  look  at  all  the  data.  Instead  the  user  can  view  a  range  of 
elements.  Thus,  when  the  option  is  selected,  the  user  is  prompted  for  the  first 
and  last  elements  in  the  range.  The  connectivities  of  these  elements  are  then 
written  to  the  screen. 

4.2.2.4  Option  4:  Optimize  Node  Numbering 

This  option  renumbers  the  nodes  so  as  to  minimize  the  bandwidth  of  the 
problem.  While  this  option  does  not  affect  the  BUCKDEL  result,  it  is 
recommended  that  it  be  utilized  as  the  problem  will  require  less  memory  and 
execute  faster.  The  user  should  keep  in  mind  that  the  node  numbering  of  the 
model  imported  via  Option  1,  will  be  altered. 

4.2.2.5  Option  5:  Return  to  Main  Menu 

Select  this  option  in  order  to  return  to  the  Main  Menu. 

4.2.3  BUCKGEN  Input  Menu  3:  Exit  BUCKGEN 

Input  Menu  3  allows  the  user  to  exit  the  BUCKGEN  program  and  save  the 
generated  data  in  the  file  specified  by  the  user  in  Input  Menu  1  (by  entering  a  1). 
The  user  can  also  use  this  option  to  quit  the  program  without  saving  the  data  (by 
entering  a  -1). 


4.3  Problem  “Dependent”  Control  Parameters  Not  Set  by  BUCKGEN 
BUCKGEN  takes  care  of  the  majority  of  BUCKDEL  input  data  file  creation  for 
the  user,  including  the  specification  of  problem  “independent”  control 
parameters.  However,  there  are  several  problem  “dependent”  control  parameters 
which  the  user  must  set  manually.  This  section  lists  these  parameters.  Reference 
is  made  to  the  input  data  file  cards  which  are  described  in  Section  5. 
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1 .  All  parameters  on  Card  5 . 

2.  The  parameter  ECCEN  on  Card  7. 

3.  The  parameters  Nil  and  N22  on  Card  8,  if  the  problem  contains  a 
delamination. 


5  BUCKDEL  Input  File  Format 

This  section  describes  the  BUCKDEL  input  file  format.  The  file  format  is 
described  as  a  series  of  cards.  In  addition  to  the  cards  described  in  this  section, 
the  input  file  may  also  contain  any  number  of  comment  cards.  These  comment 
cards  must  begin  with  ****’.  As  discussed  in  the  previous  section,  BUCKGEN 
creates  this  file  for  the  user.  Only  the  problem  “dependent”  control  parameters 
described  in  Section  4.3  must  be  manually  specified  by  the  user.  However,  the 
user  may  wish  to  change: 

•  material  properties, 

•  ply  orientations, 

•  solution  control  parameters, 

•  output  control  parameters,  etc. 

of  a  specific  problem.  To  do  so  the  user  would  need  to  either  create  a  new 
NASTRAN  bulk  data  file  for  translation  by  BUCKGEN,  or  edit  an  existing  data 
file.  If  the  latter  method  is  chosen,  the  information  contained  in  this  section  will 
be  useful. 
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Columns 


Format 


Card  1 


Quantity 


1-5 

Analysis  option  (IANALCD) 

=  1  Linear  analysis 

=  2  Linear  and  buckling  analysis 

=  3  Buckling  and  post-buckling  analysis 

15 

=  4  Nonlinear  analysis  through  both  limit  and  bifurcation  points 
=  5  Nonlinear  analysis  through  limit  points  only 

6-10 

Number  of  nodes  (NUMNP) 

15 

11-15 

Number  of  plate  elements  (NEQ) 

15 

16-20 

Number  of  beam  elements  (NEL) 

15 

21-25 

Number  of  concentrated  (nodal)  load  cards  (JSO) 

15 
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Columns _ Quantity _ _ Format 

1-5  Option  between  Pure  Newton-Raphson  and  modified  15 

Newton-Raphson  methods  (NEWTON) 

=  0  All  iterations  modified  Newton-Raphson 
&  0  first  NEWTON  iterations  with  pure  Newton-Raphson 
and  after  that  modified  Newton-Raphson  iterations 


6-10  Maximum  number  of  iterations  in  each  step  (MAXR)  15 

11-15  Maximum  number  of  steps  before  bifurcation  (MAXI)  15 

16-20  Maximum  number  of  steps  after  bifurcation  (MAXPB)  15 

21-25  Desired  (expected)  number  of  iterations  (10)  15 


26-30  Maximum  number  of  iterations  in  buckling  analysis  (ITN)  15 
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Columns 


Format 


Card  3 


Quantity 


I- 10  Given  increment  of  load  factor  at  bifurcation  point  F 1 0. 3 

(DSMITA)  (default  =  0.2) 

II- 15  Options  for  bifurcation  path  following  (KUD)  15 

=  1  follow  the  ascending  post-buckling  branch 
=  -1  follow  the  descending  post-buckling  branch 

16-25  Coefficient  in  the  formula  for  switching  post-buckling  F10.3 

branches  (COE A)  (default  =1) 


The  first  load  point  after  bifucation  ( A^)  is  determined  from  the  user  supplied  value  for 
DSMITA.  When  nonlinear  analysis  through  bifurcation  points  (IANALCD=4)  is 
performed. 


DSMITA= 


Ah  —  A, 


cr 


^ cr 

For  buckling  and  post-buckling  analysis  (IANALCD=3), 
A fr  ~  Xr 


DSMITA= 


‘'cr 


* cr 


COEA  is  used  to  determine  the  displacement  increment  ( A14,  =  14,  —  as  given  by 

AUfy  =  COEA  x  A  x  Lfeucfc  for  IANALCD=3 

AUfj  =  DSMITA  x  (14+1  -  14)  +  COEA  x  A  x  I4,udc  for  IANALCD=4 

I4,ucfcis  the  buckling  mode  and  “A”  is  determined  by  the  analysis. 
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f - : - : 

Card  4 

.  .  . 


Columns _ Quantity _ Format 


1-10 

Load  step  factor  for  the  initial  step  (DDZ1) 

F10.3 

11-20 

Tolerance  for  the  generalized  deflection  (EE1) 

F10.3 

21-30 

Tolerance  for  the  load  factor  (EE2) 

F10.3 

31-40 

Maximum  allowable  increment  of  the  load  factor  (SKLOAD)  FI 0.3 
load  increment  <  SKLOAD  x  DDZ1  at  any  load  step 

41-45 

Arc-length  method  (KARC) 

=  0  sphere  arc-length  (Crissfield) 

=  1  plane  arc-length  (Riks) 

=  2  load  control  method 

15 

The  convergence  conditions  in  the  arc-length  are  satisfied  if 


{F}r{A^-Agr‘} 

{ffK} 


<  EE1 


and 


ax.  -ax:' 


AX. 


<  EE2 


F ,  A qln  and  Attn  are  the  load  vector,  displacement  vector  increment  and  load 
vector  increment,  respectively,  at  the  n’th  load  step  and  i’th  iteration  loop. 
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Card  5 


Columns _ Quantity _ Format 


1-5 

Load  step  intervals  at  which  to  print  displacements  (IPP1) 

15 

6-10 

Node  ID  for  printing  displacement  (NPPP1) 

15 

11-15 

Degree  of  freedom  ID  for  NPPP1  (NDDD1) 

15 

16-20 

Node  ID  for  printing  displacement  (NPPP2) 

15 

21-25 

Degree  of  freedom  ID  for  NPPP2  (NDDD2) 

15 

26-30 

Node  ID  for  printing  displacement  (NPPP3) 

15 

31-35 

Degree  of  freedom  IDfor  NPPP3  (NDDD3) 

15 

36-40 

Node  ID  for  printing  energy  release  rate  (NPPP4) 

15 

The  degree  of  freedom  ID  are  1  (^-displacement),  2  (y-displacement) ,  3  (z~ 
displacement),  4  (y-rotation)  and  5  (^-rotation).  NPPP1,  NPPP2  and  NPPP3  are 
the  node  numbers  where  the  displacements  are  printed  out  in  the  output  file 
"’filename. path”.  For  example,  NPPP1=273  and  NDDD1=3,  prints  the  z- 
displacement  of  node  273.  NPPP4  is  a  node  on  the  delamination  front  where  the 
energy  release  rate  is  printed  in  the  output  file  "filename.zgm” . 
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Columns 


Format 


Card  6 


D 


Quantity _ 

1-13  Elasticity  modulus  in  fiber  direction  of  lamina  (El) 

14-26  Elasticity  modulus  in  transverse  direction  of  lamina  (E2) 

27-3 1  Poisson’s  ratio  in  plane  1-2  of  lamina  (PS  1) 

32-36  Poisson’s  ratio  in  plane  2-1  of  lamina  (PS 2) 

37-49  Shear  modulus  in  plane  1-2  of  lamina  (G12) 

50-62  Shear  modulus  in  plane  1-3  of  lamina  (G13) 

63-75  Shear  modulus  in  plane  2-3  of  lamina  (G23) 


E13.7 

E13.7 

F5.2 

F5.2 

E13.7 

E13.7 

E13.7 
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Card  7 


3 


Columns _ Quantity _ 

1-12  Axial  (extensional)  rigidity  of  beam  (E A) 

1 3-24  Bending  rigidity  of  beam  in  the  plane  perpendicular  to  plate 

(EJ) 

25-36  Bending  rigidity  of  beam  in  the  plane  parallel  to  plate  (EJ1) 

37-48  Torsional  rigidity  of  beam  (GJ) 

49-60  Transverse  shear  rigidity  of  beam  (GA) 

61-72  Feature  not  available,  set  to  0.0 

73-77  Eccentrity  (ECCEN) 

ECCEN  is  the  z-coordinate  of  the  beam  reference  axis  relative  to  the 
plane  of  the  plate  elements  (z  =  0). 


Format 

E12.5 

E12.5 

E12.5 

E12.5 

E12.5 

E12.5 

F5.2 

reference 
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^  Card  8  ) 


Columns _ Quantity _ Format 

I- 5  Number  of  plys  in  the  undelaminate  (NT)  15 

6-10  Number  of  plys  in  the  delaminate  (NT  1)  15 

II- 15  Number  of  plys  in  the  the  base  (NT2)  15 

16-20  First  undelaminate  node  on  the  delamination  front  if  the 

delamination  front  intersects  a  boundary.  If  the  delamination 


front  does  not  intersect  a  boundary,  any  node  on  the 
undelaminate  delamination  front  can  be  specified.  (Nil)  15 

21-25  Next  undelaminate  node  on  the  delamination  front  in 

the  CCW  direction  from  node  Nil  (N22)  15 

26-35  Total  laminate  thickness  (HC2)  F10.5 

36-45  z-coordinate  of  mid-surface  of  delaminate  with  respect  to  FI 0.5 

undelaminate  ( ) 

46-55  z-coordinate  of  mid-surface  of  base  with  respect  to  F 1 0.5 

undelaminate  (Z^) 

56-65  Feature  not  available,  set  to  0.0  (AK1)  F10.5 

66-75  Feature  not  available,  set  to  0.0  (AK2)  F10.5 


delaminate  (NT1  layers) 
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Columns 


Quantity 


Format 


I- 10  Ply  thickness  (T)  F10.5 

II- 20  Ply  angle  (degrees)  measured  from  the  global  x-axis  (ZT A)  F10.5 

These  ply  cards  are  written  from  the  top  of  the  laminate  to  the  bottom  (i.e.  in  the 
+z  to  -z  direction). 
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Define  the  NUMNP  nodes  and  the  displacement  and  rotational  boundary 
conditions  in  the  next  NUMNP  cards. 


Columns 

1-5 

6-10 


11-30 

31-50 

51-70 

71-75 


Cards  for  Nodal  Coordinates 


1 


Quantity 


Format 


Node  number 

15 

Code  for  displacement  boundary  conditions 
=  0  no  constraints 
=  1  x  displacement  is  constrained 
=  2  y  displacement  is  constrained 
=  3  z  displacement  is  constrained 
=  4  x  and  y  displacements  are  constrained 
=  5  y  and  z  displacements  are  constrained 
=  6  x  and  z  displacements  are  constrained 
=  7  x,y  and  z  displacements  are  constrained 

15 

^-coordinate 

E20.10 

y-coordinate 

E20.10 

z-coordinate 

E20.10 

Code  for  rotational  boundary  conditions 
=  0  no  constraints 
=  1  jc-axis  rotation  is  constrained 
=  2  y-axis  rotation  is  constrained 
=  4  Jt-axis  and  y-axis  rotations  are  constrained 

15 
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Define  the  NEQ  plate  element  connectivities  and  their  locations  in  the  next  NEQ 
cards. 


[Cards  for  P!a,e  Element^ 


Columns 


Quantity 


Format 


1-5 

Plate  element  number 

15 

6-10 

Location  of  plate  element 
=  1  delaminate  region 
=  2  base  region 
=  3  undelaminate  region 

15 

11-15 

node  1 

15 

16-20 

node  2 

15 

21-25 

node  3 

15 

For  plate  elements  on  the  delamination  front,  the  delaminate  or  base  elements  on 
one  side  of  the  front  and  undelaminate  elements  on  the  other  side,  should  be 
symmetric  about  the  delamination  front,  as  shown  in  the  figure  below.  This 
restriction  stems  from  the  energy  release  rate  calculation  (see  the  Theory  Manual 
for  more  information).  If  the  user  is  not  interested  in  the  energy  release  rates, 
this  requirement  can  be  neglected. 


delamination  front 
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Define  the  NEL  beam  element  in  the  next  NEL  cards. 


Cards 


Columns 


Quantity 


Format 


1-5 

Beam  element  number 

15 

11-15 

Node  1 

15 

16-20 

Node  2 

15 
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Define  the  JSO  concentrated  (nodal)  load  cards. 


1 


Cards  for  Concentrated  (Nodal)  Loads 


Columns 

Quantity 

Format 

1-6 

Node  ID 

16 

7-18 

Force  in  the  x-direction 

F12.5 

19-30 

Force  in  the  y-direction 

FI  2.5 

31-42 

Force  in  the  ^-direction 

F12.5 

43-54 

Moment  in  the  x-direction 

F12.5 

55-66 

Moment  in  the  y-direction 

F12.5 
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6  Output  Files 

The  output  file  names  are  based  on  the  name  of  the  input  file.  The  user  is 
prompted  for  the  input  file  name  when  BUCKDEL  is  executed.  If  the  input  file 
is  “filename,  dat”,  then  the  following  set  of  output  files  are  created. 

filename. out.  Echoes  the  main  control  variables,  nodal  coordinates,  shell  and 
beam  element  connectivities  and  material  properties  from  the  input  file.  The  file 
contains  information  about  nodes  and  elements  on  the  delamination  front.  The 
linear  solution  due  to  reference  load,  buckling  mode  shape  and  buckling  load  are 
also  available. 

filename, femesh:  The  data  in  this  file  can  be  used  to  plot  the  finite  element  mesh 
of  the  model  using  GNUPlot. 

filename. buck:  The  data  in  this  file  can  be  used  to  plot  the  fundamental  buckling 
mode  shape  using  GNUPlot. 

filename. zz\  The  energy  release  rate  at  each  node  on  the  delamination  front  at 
each  load  step. 

filename. 7.pvn:  The  load  factor  (1st  column),  maximum  energy  release  rate  (2nd 
column)  and  the  energy  release  rate  at  node  NPPP4  (3rd  column),  at  each  load 
step. 

filename. path:  The  load  factor  (1st  column),  displacements  at  nodes  NPPP1 
(2nd  column),  NPPP2  (3rd  column),  NPPP3  (4th  column)  at  each  load  step. 

/i/gname.diff:  Load  factor  and  norm  of  the  displacements  along  z-direction  of 
nodes  belonging  to  the  base  and  undelaminate  regions  at  each  load  step. 

zmode.xx:  z-displacements  for  plotting  with  GNUPlot  at  load  step  intervals 
specified  by  IPP1  (card  5,  field  2).  The  extension  “xx”  is  determined  by  the 
interval,  IPP1.  For  example,  if  IPP1  is  4,  then  a  sequence  of  files,  “zmode.04”, 

“zmode.08”,  “zmode.12”, . ,  etc.,  are  created  containing  the  z-displacements 

at  the  specified  interval. 

Files  “filename. zg”,  “filename. zgm”,  “filename. diff’  and  “filename.pa.th”  are 
empty  when  the  linear  analysis  options  (IANALCD  =  1  or  2)  are  specified. 

To  view  the  plot  files  “filename. buck”  or  “zmode.xx”,  execute  GNUPlot  and 
type  the  following  at  the  prompt. 

>  set  parametric  off 

>  splot  “filename. buck”  w  1 
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To  view  the  finite  element  mesh  in  “filename. femesh”,  type  the  following  at  the 
prompt 

>  plot  “filename. femesh”  w  1 

GNUPlot  can  be  used  to  create  plots  from  the  data  files  “filename. path”  or 
filename. zgm” .  For  example,  if  it  is  desired  to  plot  the  displacement  at  node 
NPPP2  (column  3)  vs.  load  factor  (column  1)  from  data  file  “ filename: path”, 
type  the  following  at  the  prompt, 

>  plot  “filename. path”  u  1:3  w  1 

Similarly,  to  plot  the  norm  of  the  ^-displacements  of  the  base  and  undelaminate 
regions  vs.  load  factor,  type 

>  plot  “filename. diff  ’  u  1:2  w  1 


7  Example  Problems 

This  section  will  illustrate  some  of  the  capabilities  of  BUCKDEL  by  solving 
several  example  problems.  The  input  and  output  files  for  each  example  are 
distributed  with  the  software  so  that  the  user  can  utilize  them  in  learning  how  to 
run  BUCKDEL. 

7. 1  Simply  Supported  Laminated  Plate  Under  Biaxial  Compression 
A  simply  supported  graphite-epoxy  laminated  plate  under  biaxial  compression  is 
shown  in  Figure  7-1.  The  loads  per  unit  length  in  the  x  and  y  directions  are  ANX 
and  ANy ,  respectively,  with  X  being  the  load  factor.  The  lamina  properties  are: 
Ei=18.5e+6  psi,  E2=1.89e+6  psi,  Gi2=0.93e+6  psi,  and  v12=0.3.  The  laminate 
has  16  plies  with  a  ply  thickness  of  0.005  inches  (total  laminate  thickness  =  0.08 
inches).  The  buckling  load  for  several  stacking  sequences  and  values  of  Ny/Nx 
(Nx=1.0  psi)  was  computed  with  BUCKDEL.  Due  to  symmetry,  only  one 
quarter  of  the  plate  was  modeled  with  finite  elements  (see  Figure  7-2).  The 
results  are  shown  in  Table  7-1  where  they  are  compared  to  those  contained  in 
[1].  In  general,  the  agreement  is  very  good.  The  slight  differences  are  due  to  the 
fact  that  BUCKDEL  accounts  for  shear  deformation  while  the  results  from  [1]  do 
not.  Figure  7-3  is  a  plot  of  the  fundamental  buckling  mode  (Ny/Nx  =  2.0  case) 
produced  by  using  GNUPlot  to  post-process  the  “*.buck”  file  created  by 
BUCKDEL. 
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Table  7-1:  Buckling  loads  for  several  stacking  sequences  and  load  ratios 


Ny/Nx 

Stacking  Sequence 

A  cr  (BUCKDEL) 

A  cr  (Ref. 

[1]) 

0.2 

(45,-452,90,452,90,-45)s 

134.82 

137.10 

1.0 

(90,45, 904,-45,90)s 

56.18 

61.99 

2.0 

(905,45, -45, 90)s 

35.68 

36.84 

XNy 


ANy 

Figure  7-1:  Simply  supported  graphite-epoxy  laminated  plate  under  biaxial  compression 
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’exOla.buck’  - 


Figure  7-3:  Fundamental  buckling  mode  (Ny/Nx  =  2.0  case) 


7.2  Simply  Supported  T  .aminated  Plate  with  Delamination 
The  simply  supported  graphite-epoxy  laminated  plate  in  Section  7.1  is 
reanalyzed,  except  that  in  this  section  a  delamination  is  introduced.  The  lamina 
properties,  geometry,  etc.  are  the  same  as  in  Section  7.1.  The  delamination  is 
centrally  located  in  the  plate.  The  load  ratio  is  Ny/Nx  =1.0  and  the  ply  stacking 
sequence  is  (90,  45,  904,  -45,  90)s.  The  analysis  was  performed  with  the 
delamination  located  at  the  mid-surface  of  the  laminate  (i.e.  between  the  90 
plies)  and  between  the  2nd  and  3rd  plies  (45  and  90  plies)  from  the  top  face  of 
the  plate.  The  buckling  load  with  the  delamination  at  the  mid-surface  is  40.46 
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and  the  corresponding  buckling  mode  can  be  seen  in  Figure  7-4,  which  was 
produced  by  using  GNUPlot  to  post-process  the  “*.buck”  file.  When  the 
delamination  is  located  between  the  2nd  and  3rd  plies,  the  buckling  load  is  15.09 
and  the  buckling  mode  is  plotted  in  Figure  7-5. 

To  demonstrate  the  difference  between  a  “global”  and  “local”  buckling  mode 
(see  the  Theory  Manual),  it  is  useful  to  solve  the  same  problem,  except  with  the 
delamination  replaced  by  a  hole  having  the  same  dimensions  as  the 
delamination.  If  this  is  done,  the  buckling  load  is  found  to  be  23.80  and  the 
buckling  mode  is  plotted  in  Figure  7-6.  By  comparing  Figures  7-4,  7-5  and  7-6 
one  can  observe  that  the  buckling  mode  in  Figure  7-5  is  “local”  while  those  in 
Figures  7-4  and  7-6  are  “global”.  The  importance  in  making  this  distinction  is 
that  if  the  user  wants  to  investigate  the  effect  of  a  delamination  on  the  global 
buckling  load,  it  is  often  necessary  to  perform  a  nonlinear  post-buckling  analysis. 
A  linear  buckling  (eigenvalue)  analysis  of  a  structure  with  a  delamination  will 
often  yield  the  “local”  buckling  load.  The  structure  will  still  have  a  great  deal  of 
load  carrying  capability  left  after  “local”  buckling  occurs.  Thus,  a  design  based 
on  the  “local”  buckling  load  would  be  overly  conservative. 


' ex02a.buck'  - 


Figure  7-4:  Fundamental  buckling  mode,  delamination  at  mid-surface 
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'ex02b.buck' 


4.! 


Figure  7-5:  Fundamental  buckling  mode,  delamination  between  2nd  and  3rd  plies 


'ex02c.buck'  - 


Figure  7-6:  Fundamental  buckling  mode,  plate  with  a  hole 


8  References 

[1]  J.L.  Walsh  and  R.T.  Haftka;  Stacking-Sequence  Optimization  for  Buckling  of 
Laminated  Plates  by  Integer  Programming,  AIAA  Journal,  Vol.  30,  No.  3,  pp. 
814-819,  March  1992. 


26 


Figure  1:  Multi-domain  model  of  a  laminated  plate  containing  a  delamination 

1  INTRODUCTION 

Laminated  composites  are  gaining  importance  in  aircraft  structural  applications  due 
to  their  very  high  strength-to-weight  ratios.  These  applications  are,  however,  still  very 
limited  due  in  part  to  the  lack  of  design  tools  capable  of  modeling  the  failure  mechanisms 
involved  in  the  laminated  structures.  BUCKDEL  is  a  design  tool  which  addresses  this 
need.  The  BUCKDEL  software  performs  geometrically  nonlinear  analysis  of  stiffened 
laminated  composite  shells  with  and  without  delaminations. 

BUCKDEL  allows  the  user  to  perform: 

•  a  linear  static  solution; 

•  a  linear  buckling  (eigenvalue)  analysis;  and 

•  a  nonlinear  post-buckling  analysis  through  both  limit  and  bifurcation  points. 

There  are  two  finite  elements  implemented  in  BUCKDEL;  a  quasi-conforming  triangular 
laminated  shell  element  based  on  a  refined  first  order  shear  theory  with  seven  degrees  of 
freedom  per  node  and  a  fiber-reinforced  beam  element  with  seven  degrees  of  freedom  per 
node.  The  behavior  of  all  materials  is  assumed  to  be  linearly  elastic. 

The  multi-domain  modeling  method  is  used  for  the  analysis  of  laminated  plates /shells 
containing  delaminations.  In  this  model,  the  delaminate,  base  and  the  undelaminate 
are  modeled  as  3  distinct  plates  (see  Figure  1).  On  the  delamination  front,  Mindlin’s 
deformation  assumption  is  applied  to  obtain  the  relationship  between  the  displacements 
in  the  delaminate,  base  and  undelaminate. 

The  remainder  of  this  manual  provides  descriptions  of  the  theoretical  concepts  asso¬ 
ciated  with  the  BUCKDEL  software.  This  begins  with  some  basic  concepts  in  damage 
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mechanics  as  it  applies  to  laminated  composite  materials.  Thereafter,  a  detailed  descrip¬ 
tion  of  the  finite  elements  and  solution  schemes  implemented  in  BUCKDEL  are  given. 

2  DELAMINATIONS  IN  LAMINATED  COMPOS¬ 
ITE  MATERIALS 

It  is  known  that  delaminations  are  the  most  frequent  causes  of  failure  of  laminated 
structures,  particularly  under  compressive  loads.  The  presence  of  delaminations  leads  to  a 
reduction  in  the  overall  buckling  strength  of  the  structure.  In  addition,  the  delaminations 
tend  to  grow  rapidly  under  post-buckling  loads,  in  particular,  causing  a  further  reduction 
in  structural  strength,  and  finally,  leading  to  fatal  structural  failure. 

Fiber  reinforced  Polymer  Matrix  Composites  can,  in  general,  be  modeled  as  trans¬ 
versely  isotropic  linear  elastic  solids.  Thus,  each  lamina  in  a  laminated  composite  struc¬ 
ture  can  be  modeled  as  a  transversely  isotropic  solid.  However,  when  a  laminate  with 
an  arbitrary  (but  perhaps  optimized)  stacking  sequence  of  angle-ply  lay-up  is  consid¬ 
ered,  additional  couplings  in  the  constitutive  matrix  result.  In  aircraft  structures  made 
up  of  PMC’s,  one  can  tailor  the  strength,  stiffness,  aeroelastic  and  acoustic  response  of 
the  structure  by  optimizing  the  constituent  properties  as  well  as  the  stacking  sequence. 
However,  a  stacking  sequence  that  optimizes  the  stiffness  and  dynamic  response  may  not 
necessarily  be  optimal  from  the  point  of  view  of  controlling  the  deleterious  interlaminar 
stresses  which  may  induce  local  delaminations. 

An  obvious  consequence  of  delamination  on  the  bending  behavior  of  a  beam  or  plate 
is  the  complete  absence  of  shear  transfer  from  outer  plies  -  carrying  larger  normal  stress 
-  to  inner  plies,  except  through  friction  when  such  tangential  forces  exist.  However,  more 
critical  than  the  lack  of  the  shear  transfer,  the  composite  laminate  with  delaminated 
areas  may  have  local  buckling  phenomena  in  the  presence  of  local  compression  loads 
when  they  exceed  certain  critical  buckling  loads  at  the  ends  of  the  delamination  sites. 
Depending  on  the  area  and  location  of  the  delamination,  the  buckling  may  have  different 
modes  (Figure  2). 

The  global  behavior  of  primary  interest  when  designing  aircraft  structures  made  of 
PMC’s  is  the  buckling  and  post-buckling  response  of  stiffened  laminated  composite  plate 
and  shell  structures  that  contain  reinforced  or  unreinforced  cut-outs.  Delaminations 
and  local  buckling  may  accompany  and  concurrently  influence  global  buckling.  The 
accurate  prediction  of  the  onset  of  local  buckling  due  to  delaminations  and  its  effect  on 
the  global  buckling  and  post-buckling  response  of  a  structure  is  essential  in  establishing 
the  integrity  of  that  structure.  BUCKDEL  provides  a  means  of  modeling  the  strength 
reduction  due  to  the  presence  of  the  delaminations  and  the  mechanisms  of  delamination 
growth.  BUCKDEL  can  be  used  for  modeling  delaminated  stiffened  composite  structures; 
for  an  accurate  prediction  of  their  geometrically  nonlinear  behavior  under  normal  and 
post-buckling  loading  conditions;  and  for  delamination  growth  prediction  by  using  the 
pointwise  energy  release  rate  at  the  delamination  front. 
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Figure  2:  Buckling  modes  for  a  delaminated  composite 
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3  MULTI-DOMAIN  MODELING  OF  DELAMINA¬ 
TIONS 


BUCKDEL  can  be  used  to  analyze  a  general  stiffened  laminated  composite  plate/shell 
with  a  single  delamination  of  an  arbitrary  shape  and  location,  subjected  to  arbitrary 
compressive  loads  (Figure  3).  The  structure  is  modeled  using  the  multi-domain  model 
wherein  the  delaminated  plate/shell  is  assumed  to  be  assembled  with  three  distinct  plates 
or  shells  -  (1)  Laminate:  nondelaminated  zone  f (2)  Delaminate:  thinner  side  of  the 
delaminated  zone  and  (3)  Base:  thicker  side  of  the  delaminated  zone  The  three 
shells,  =  1, 2, 3  respectively,  have  midsurface  areas  thicknesses  boundaries 

and  midsurface  boundaries  dA^\  The  delamination  edge  is  denoted  by  T.  The 
assumptions  of  Reissner-Mindlin  theory  of  plate  bending  are  used  for  modeling  each 
plate/shell  and  the  joint  between  them.  Thus,  for  each  plate/shell,  the  3-dimensional 
displacement  field  (U  =  {U\  U2  U$})  can  be  expressed  in  terms  of  the  corresponding 
midsurface  displacement  (u  =  {uj  u2  U3})  and  rotation  (6  =  {9i  02  0})  fields  as, 

U«(a:a,X3)  =  u (i)(xa)  -  (1) 

where  xW  (a  =  1,2)  are  the  inplane  curvilinear  shell  coordinates  and  is  the  thickness 
coordinate  for  the  ith  (i  =  1,2,3)  shell  (Figure  3).  The  structural  continuity  at  the 
delamination  front  F  is  maintained  by  assuming  the  deformation  to  be  unique  at  the 
junction  of  the  three  shells  i.e.  =  U2)  =  U3)  on  F  in  accordance  with  the  Reissner- 
Mindlin  law  of  flexure  [equation  (1)].  In  other  words,  at  the  delamination  edge ,  the 
mid-surface  degrees  of  freedom  of  the  delaminate  and  the  base  shells  are  assumed  to  be 
related  to  those  of  the  nondelaminated  shell  by, 


=  42) 

=  43)  ^ 

=  *£2) 

=  *i3) 

(2) 

«s? 

= 

+  h®0W  J 

at  r 

where  is  the  distance  of  the  midsurface  of  the  ith  shell  from  the  laminate  midsurface 
(Figure  3).  It  can  be  noted  that  the  above  continuity  conditions  at  the  delamination 
edge  can  be  modified  appropriately  when  using  any  other  alternative  shell  theory  (e.g. 
higher  order  shear  deformable  theory)  or  by  choosing  appropriate  heuristic  multi-point 
constraints  based  on  experience. 

Similarly,  the  beam  (stiffener)  degrees  of  freedom  are  related  to  the  shell  degrees  of 
freedom  such  that  the  transverse  variations  of  deformation  across  the  shell  and  beam 
section  are  consistent  with  Reissner-Mindlin  theory: 


u3  —  U3 

oba  =  esQ 

ua  =  <4  +  e^a 


(3) 


where  superscripts  6  and  s  represent  beam  and  shell  degrees  of  freedom  respectively,  and 
e  is  the  eccentricity  of  the  stiffener’s  neutral  axis  with  reference  to  the  neutral  surface  of 
the  shell. 
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(a)  Del  animation  in  an  arbitrary  plate 


K - L - w 


(b)  Clamped  square  plate  with  central  elliptic 
delamination 


Delaminate  ft® 


1 


Assumed  deformation  at  Plate-Junction 

(Ud)  =  U<2>  -  U(3))r 


Figure  3:  Multi-domain  model 


If  4>a  (a  =  1,2)  are  the  midsurface-rotations  characterizing  the  transformation  of  the 
normal  to  the  undeformed  mid-surface  to  that  of  the  deformed  mid-surface,  we  have, 


<j)a  —  U3,a  ”1“  bapUp 


(4) 


where  b  is  the  curvature  tensor  of  the  shell’s  mid-surface. 

In  the  present  formulation,  the  engineering  transverse  shear  strains  7^3  are  considered 
as  independent  mid-surface  degrees  of  freedom,  assuming  that  they  do  not  vary  over 
the  shell  thickness.  The  engineering  transverse  shear  strains  are  related  to  the  section 
rotations  9a  and  the  mid-surface  rotations  <f>a  in  accordance  with  the  Reissner-Mindlin 
theory  of  plate  flexure,  as, 

7a3  =  <f>a-0a  (5) 

For  each  shell,  the  inplane  strain  components  (en,  e22,  £12)  and  the  transverse  shear 
strain  components  (ei3,  €23)  of  the  3-dimensional  Green-Lagrange  strain  tensor  (including 
large  deformations)  are  given  by, 

=  j  [Ua,f  +  C'5,„  +  04,,  +  K;,US)  CAg  +  b,„U0)  -  baeU3f 

42  =  (6) 

respectively  where,  again,  or,  /?  =  1,2  and  i  =  1,2,3.  Using  the  3-dimensional  displace¬ 
ment  field  description  [equation  (1)]  in  equations  (6),  and  neglecting  the  nonlinear  terms 
associated  with  the  section  rotations,  the  inplane  components  of  the  Green-Lagrange 
strain  tensor  can  be  simplified  as, 


where,  s^p  and  ufp  are  the  linear  and  nonlinear  components  of  the  membrane  strains; 

and  /c^g  and  Xap  are  the  flexural  strain  components  related  to  mid-surface  rotations  <fa 
and  transverse  shear  strains  70,3  respectively,  given  by, 

£a0  =  \  “  [KpU3]^ 

ual  =  \  [(«3,er  +  hf}Up)  (u3<p  +  6/3aUa)](0 

Ka0  ~  ~2~  [(U3,a  +  bapUp)  ,p  +  («3,/3  +  bp aUa)  ,a 

Xal  =  \  [7a3,/3  +  7/33, a]W  (8) 


Equation  (7)  can  be  written  in  a  compact  form  to  express  the  3-dimensional  Green- 
Lagrange  strain  components  eh)  =  {eu  e22  £12  tx3  £23}^  in  terms  of  the  2-dimensional 
strain  components  eh)  =  {en  g22  e12}h);  i/h)  =  { t/22  i72}h);  k h)  =  {kh  k22  ^12}^); 
=  {xn  X22  Xi2}W  and  7W  =  {7x3  723}h)  as, 


.(*)  _  |(£  +  ^)  +  ^3(k  +  x)\W 

i  h 


I 


=  B  S 


(9) 
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where  S  =  {ui  u2  u$  u3} i  u3i 2  713  723}  is  the  vector  of  independent  degrees  of  freedom 
and  B  is  the  strain-displacement  matrix. 

Assuming  orthotropic  material  behavior  for  each  lamina,  the  inplane  stresses  = 
{an  (722  012}^  and  the  transverse  shear  stresses  =  {r13  r23}^  are  related  to  the 
corresponding  strain  components  as, 

En  Eie  0  0 

E12  E22  E26  0  0 

Eie  E26  Eqq  0  0 

0  0  0  £44  -E45 

0  0  0  £45  £^55 

where  the  material  constitutive  terms  E^  are  functions  of  the  thickness  coordinate  of 

each  shell  x$K  Generally,  for  a  laminate  with  orthotropic  layers,  are  assumed  to  be 
piecewise  constants  over  the  laminate  thickness. 

Now,  for  each  shell,  the  strain  energy  density  (per  unit  volume)  can  be  calculated 
using  equations  (9)  and  (10)  as, 

W(,)  =  [tr  •  ((e  +  u)  +  x3  (k  +  x))]W  +  [r  ■  t]W  (11) 

As  we  can  observe  from  equations  (7)  through  (10),  the  strain  field  is  a  linear  function  in 
X3  while  the  stress  field  is  a  piecewise  linear  function  in  x3.  Hence,  it  is  often  more  useful 
to  integrate  the  strain  energy  density  explicitly  in  the  thickness  direction  and  define 
strain  energy  density  per  unit  area  as, 

=  [  W®dx. 3  =  [N  •  (e  +  v)  +  M  •  (k  +  x)  +  Q  ■  t](0  (12) 

Jti 

where  N  =  { Nn  N22  A12}  are  the  inplane  stress  resultants;  M  =  {Mu  M22  M12}  are  the 
bending  moment  resultants;  and  Q  =  {Q 13  Q23}  are  the  transverse  shear  stress  resultants 
given,  for  each  shell,  by, 

'  N  )  (,)  [A  B  0  ] (i)  f  (e  +  v)  ]  (i) 

<  M  >  =  BDO  {(k  +  x)’  (13) 

.QJ  Lo  0  Gj  It, 

and 

(AW;  Bia\  Dkl){i)  =  [  £$(*3)  (l;  x3;  xj)  dx 3 

Jti 

Gmn  =  l  SmSnE§{x3)dx3 

as  k,l  =  1,2,3  correspond  to  i,j  =  1,2,6  and  m,n  =  1,2  correspond  to  i,j  =  4,5; 
and  Si  and  S2  are  the  shear  correction  factors  in  the  transverse  planes  1—3  and  2  —  3 
respectively.  Integrating  [equation  (13)]  over  the  shell  domains  we  get  the  total 
strain  energy  stored  in  the  ith  shell  as, 

WM  =  \J{i)W®dA  =  i(N  •(£  +  !/)  +  M-(K  +  x)  +  Q-7)(t)^A  (14) 
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If  f(0  =  {/,  /2  /3}(0  are  the  body  force  components  acting  on  the  ith  shell;  and 
F  =  {F1  F2  F3}«  are  the  traction  forces  acting  on  the  boundaries  of  the  ith  shell  -  all 
defined  in  the  curvilinear  shell  coordinate  system  -  the  potential  of  the  external  forces  on 
the  ith  shell  is  given  by, 

=  f  f  •  U  <m+  f  F-U  dA  (15) 

Assuming  that  the  body  force  components  and  the  transverse  traction  force  component 
do  not  vary  over  shell  thickness  and  the  inplane  traction  force  components  vary  linearly 
over  the  thickness,  we  get, 

P(i>  =  L  (f-“)W"  +  L  (F'“  +  M.0)(,)dr  (16) 

where, 


and  =  {Mi  M2  0}^  represent  the  traction  moments  acting  on  dA^\ 

Then,  using  equations  (14)  and  (15),  the  total  potential  energy  for  the  assembly  of 
shells  is  given  by, 

n  =  J2  (^(i)  -  p{i))  (1?) 

1=1 

Applying  the  minimum  total  potential  energy  principle  (<HI  =  0)  ,  we  obtain  the  follow¬ 
ing  internal  equilibrium  equations  for  each  shell, 

(  (Nap  -f-  barjMrip)  ip  barjNrj0(f>0  4"  fa  0  \ 

I  Ma0tp  +  Qa 3  =  0  I  (IS) 

V  (Qa3  4"  Nap  (u3,0  4"  b pypJejf) )  , a  4-/3  0  /  ^(t) 


and  the  following  boundary  conditions  on  the  external  boundary  of  each  shell, 


I  ^3, a  4"  barjUrj  'J a3 
V  «3 


0  or 
0  or 
0  or 


Nap-Tip 

Map.np 

(Qa 3  4-  Nap  (uz,p  4-  bpvuv ))  .na 


Fa  \ 

Ma 

A 

)  PAW 


(19) 


where  n,  is  the  ith  component  of  the  unit  vector  n  normal  to  delamination  front. 
The  displacement  is  continuous  at  the  delamination  edge  and  therefore,  on  F  we  have: 


(<su(1)  =  <su(2)  =  £u(3))r 


(20) 


Hence,  for  equilibrium,  the  following  conditions  have  to  be  satisfied  at  any  point  on  the 
delamination  edge: 


{JFn(N)  =  0  ;  ^n(M')  =  0  ;  Tn(T)  =  0}r  (21) 

where  Fn(*)  =  (*)(1)  -  (*)(2)  -  (*)(3)  with  (*)(i)  corresponding  to  the  value  of 

(*)  at  the  specified  point  on  the  delamination  boundary  of  the  ith  shell;  (M7)^  = 


34 


MW  +  h® N(i);  and  Ta  =  Qa 3  +  Na0  («3)/ 3  +  b^u^)  is  the  effective  shear  force.  It  can 
be  noted  that  these  conditions  are  normally  satisfied  in  a  finite  element  model.  However, 
conflict  may  arise  if  the  delamination  edge  touches  the  structural  boundary,  particularly 
the  boundary  with  specified  displacements  where  the  stresses  are  expected  to  develop. 
Since  the  delamination  generally  acts  as  an  imperfection,  it  is  often  possible  that  the  base 
can  come  in  contact  with  the  delaminate  after  global  buckling  sets  in,  and  then  the  above 
equilibrium  equations  (21)  may  not  be  valid.  Thus  the  model  can  be  effectively  used  for 
examining  the  growth  of  embedded  delaminations  under  post-buckling  conditions  as  long 
as  the  delaminated  plies  do  not  come  in  contact  with  each  other. 

4  FINITE  ELEMENT  FORMULATIONS 

A  2-noded  curved  beam  element  (BEAM2)  and  a  3-noded  shell  element  (SHELL3)  are 
used  by  BUCKDEL  for  modeling  stiffened  shells.  The  transverse  shear  deformation  is 
introduced  explicitly  in  accordance  with  the  Reissner-Mindlin  plate  theory.  Reduced 
integration  is  used  for  the  constrained  membrane  strain  energy  component  to  eliminate 
locking  from  the  elements.  The  beam  degrees  of  freedom  are  related  to  the  corresponding 
shell  degrees  of  freedom  in  accordance  with  the  Reissner-Mindlin  theory  of  flexure. 

4.1  BEAM2:  2-noded  curved  beam  element 

The  element  incorporates  the  following  independent  degrees  of  freedom:  u  -  inplane 
displacement;  w  -  transverse  displacement;  w,x  -  slope  of  midsurface  deflection  and  7  - 
the  transverse  shear  strain.  Thus  the  element  requires  C^-continuous  fields  for  u  and  7, 
and  (^-continuous  field  for  w.  Linear  Lagrangian  polynomials  are  used  for  geometric 
description  x  and  for  the  inplane  displacement  component  u: 


x  =  hi  +  h2(  (22) 

u  =  Qi  +  a2f  (23) 

and  cubic  Hermitian  polynomials  are  used  for  the  transverse  deflection  w: 

w  =  /?i  +  /?2£  +  PsC  +  (24) 

where  £  is  the  natural  coordinate;  and  hi,  ot  and  Q%  are  functions  of  nodal  values  Xi,  Ui 
and  Wi  respectively. 

4.2  SHELL3:  3-noded  triangular  curved  shell  element 

The  element  is  described  in  the  curvilinear  coordinate  system  x  —  y  and  the  area  coordi¬ 
nates  are  used  for  field-description.  Accordingly  we  have, 

=  (25) 

{=  1 
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(26) 


Inverting  the  above  relationship,  we  get, 

Li  —  2^  (ot'l2  +  +  Ci3) 

where 

an  =  yj  yk\  cii2  =  Xk  d{  3  =  Xjyk  x^yj 

A  =  ^  (X2V3  -  +  *3yi  -  xm  +  xm  -  x2yi ) 

and  j  =  2,3,1;  k  =  3, 1,2  as  i  =  1,2, 3. 

The  inplane  displacements  and  the  transverse  shear  strains  need  to  satisfy  C°  -  con¬ 
tinuity  while  the  transverse  deflection  need  to  satisfy  (^-continuity  in  the  present  for¬ 
mulation.  The  independent  field  variables  u,  v,  w,  jxz  and  ~jyz  are  expressed  in  terms  of 
the  nodal  degrees  of  freedom  it*,  u,,  W{,  gX{  =  {-w,y ),-,  gyi  =  (w,x ),-,  ^XZi  and  jyZi  as, 


3 

{u  V  7X2  ''iyz }  —  'y  ]  Li{u  V  7x2  7y*}t 
*=1 
3 

W  =  J2(NiiWi  +  N2iQXi  +  N3igyi)  (27) 

2=1 

where, 

Nu  =  Li  +  LfLj  +  L]Lk  -  LiL)  -  L{L2k 

N2i  =  Ufci  Lj  -j-  -LiLjL^j  4*  dji  (^L{ L^  -f-  —LiLjL^j 

LT3i  =  cih2  (jL^Lj  +  —LiLjLk'j  +  aj2  ^LjLk  +  —LiLjLi^J  (28) 

are  the  cubic  polynomials  for  the  transverse  deflection. 

The  quasi-conforming  shell  element  (Atluri,  Huang  and  Shenoy  [1])  is  developed  from 
the  Hu-Washizu  variational  principle  as  follows.  If  the  inter-element  C°  continuity  is 
satisfied  a  priori,  the  Hu-Washizu  principle  can  be  written  simply  as  the  sum  of  the 
respective  integrals  for  each  of  the  finite  elements  (Atluri  [2]),  as: 


{A(e)  +  <r:  ^(Vu  +  uv  +  Vu •  uV)  —  e  -f°-ujdy 

—  f  t°  •  ucL4  +  f  n  •  a  •  (u  —  u°)dA\  =  0  (29) 

^  Sam  •*  Sum  J 


where  a  and  e  are  respectively  the  second  Piola-Kirchhoff  stress  tensor  and  Green  strain 
tensor,  Sum  and  Scm  are  the  segments  of  the  boundary  of  element  m  where  displace¬ 
ments  and  tractions  are  prescribed  respectively.  If  u  =  u°  on  SUm,  and  a  is  considered 
as  a  Lagrange  multiplier,  Eq.  (29)  may  be  cast  into  the  following  conditional  two-field 
variational  problem 


SIl(f  {A{e)-i°-u}dV-  f  t°-udA)=0 

m  \vJ2m  ^  Sam  ' 


(30) 


-(Vu  +  uV  +  Vu  •  uV)  —  e  =  0 

/j 


(m  =  1,2,  ....N) 


(31) 
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for  all  element  strain  and  displacement  fields,  where  N  is  the  number  of  elements.  The 
constraint  condition  in  Eq.  (31)  may  now  be  satisfied  in  a  weak  form  if  the  independently 
assumed  strain  fields  e  are  expressed  in  terms  of  the  element  displacement  fields  through 
the  element-level  “quasi-conforming”  condition: 


/ 

J  £ln 


5a  : 


L2 


(Vu  +  uV  +  Vu  •  uV)  —  e 


dV 


0 

(m  =  1,2,  ....iV), 


(32) 


8a  is  an  independent  test-function  for  each  element.  When  the  element  strain  fields  e  are 
expressed  in  terms  of  element  displacement  fields  through  Eq.  (32),  the  two-field  problem 
reduces  to  a  one-field  variational  problem: 


8J2  (  f  (A(e(u))  -  f°  •  u  }dV-  f  t°  •  u  dA)  =  0 

jjl  ^  Sam  / 

(33) 

U  =  U°  on  Sum 


The  strain  energy  of  the  laminated  composite  shell  element  can  be  expressed  in  the 
form 

[  A(e)dV  =  -  [  (Na/3caPo  +  M“V/3  +  Qa 7a/?)  dA ,  (34) 

J £lrn  "  ''•Am. 

where  Am  is  the  area  of  the  element. 

Since  8a  in  Eq.  (32)  is  a  continuously  differentiable,  but  otherwise  arbitrary,  test 
function  that  varies  in  the  three  co-ordinates  xz  of  the  shell  space,  we  assume  a  specific 
form  for  such  a  test  function,  denoted  by  e*,  possessing  a  linear  variation  in  the  thickness 
direction  x3.  We  thus  have 


e*  =  (c*«0  +  +  e*“3aaa3. 

Using  Eq.  (35)  in  Eq.  (32)  and  integrating  through  x3,  we  obtain 
*0/3  f  1 


h  e*al 3  [( ua;/3  +  up.a )  -  2baPw  +  (w,a  +  by„){wtp  +  6£tv)  -  ea/3o]  | 

+  T2L  ^  ^W'°  +  +  ^  +  bPu^’a  ~  7a3;/3  ~  7/33;q]  ~  Ka/3} 

+  h  f  e*a3  (wta  -f  5£iv  -  0Q  -  7a3)  dA  =  0. 

J  Am 

Furthermore,  tap0  and  Kap  can  be  divided  into  two  parts,  respectively  as 

e  *  =  M  +e{2l 


(35) 


dA 

dA 

(36) 


(37) 


(38) 
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where  are  the  components  of  the  linear  strain  tensor,  k'q(3  is  the  total  curvature 
change  of  the  middle-surface, 

+  u/3,a)  -  bapw  (39) 


K'a/3  =  -^[(w,a  +  6>M);j5  +  (u;)jS  +  6^);a],  (40) 

(2) 

and  cap0  and  k"^  denote  the  non-linear  terms  of  the  in-plane  strain  tensor,  and  curvature 
change  due  to  transverse  shear  alone,  respectively: 


^apo  2 

K'ap  =  2  ^a3;/?  "*■  T/?3;a) 


(41) 

(42) 


If  Eq.  (42)  are  satisfied  a  priori,  sufficient  conditions  for  the  satisfaction  of  Eq.  (36)  follow 


+ 


L  e*ap  + up-'a  -  2ba*w)  -  dA 

Sa  e*a/?  Ku,,“  +  bau^(W’P  +  bPU»)\  -  dA  =  0 


(43) 


+  Kun);P  +  iw,P  +  ^puu);a 


dA  =  0 


(44) 


/  +  =  0,  (45) 

J  Am 

where  the  notation  e*a3  is  replaced  by  <f)*a.  When  Eq.  (45)  is  satisfied  for  a  complete 
family  of  test  functions  or  the  areas  Am  are  infinitesimal,  we  have 

(u^,a  -f-  <j>a  —  0. 

Thus  Eq.  (43)  reduces  approximately  to 

JA  e  1 2^Ua'A  uP<a  ~  2bapw)  —  ei(90|  dA  =  0;  (46) 

when  Eq.  (45)  is  satisfied  for  a  cut-off  family  of  test  functions  (including  a  constant  one) 
and  the  areas  Am  are  sufficiently  small. 

The  Eqs.  (46),  (44)  and  (45),  as  the  “quasi-conforming”  conditions,  are  used  to  de¬ 
termine  the  undetermined  parameters  of  e^j0,  K'aP,  4>a  in  terms  of  the  displacement  pa¬ 
rameters,  which  are  subsequently  substituted  into  Eq.  (33). 

In  the  quasi-conforming  element,  the  linear  part  of  strains  e^j0,  the  total  curvature 
change  of  the  middle  surface  k'^  and  rotations  </>Q  in  Eq.  (33)  are  discretized  directly 
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and  the  undetermined  parameters  in  them  are  expressed  in  terms  of  the  displacement 
field,  viz.  the  generalized  displacements  at  the  nodes.  In  addition,  the  transverse  shear 
strains  7^3  in  Eq.  (33)  are  discretized  by  means  of  their  values  at  the  nodes.  Based  on  a 
Lagrangean  description,  using  Eq.  (33),  the  formulation  of  the  incremental  equilibrium 
equations  in  terms  of  the  generalized  displacements  is  now  presented. 

Let  (A;,  u i)  denote  a  point  (i)  on  the  equilibrium  path,  (A;+i,  u,-+i)  is  a  required  point 

(i  +  1) 


Ui+i  =  U;  +  All; 


A;+i  —  A;  +  A  A; 


(47) 


where  A  is  the  load  factor. 


Substituting  Eqs.(37,  38  )  into  the  constitutive  relation  for  the  composite  laminate, 


we  have 


/  n  \  _  r  nw  )  r  n(2)  \ 

\  M  J  \  MW  J  +  \  M<2>  J 


(48) 


where 

(  J  N(1)  \  /  N<2>  \  \  _  [  A  B 
^  \  M«  J  ’  \  MW  J  J  “  [  B  D 

and  N  and  M  are  the  stress  resultants  and  couples  and  the  terms  A,  B,  and  D  involve 
the  laminate  material  constants. 

Let  77  and  A77  denote  the  generalized  nodal  displacements  corresponding  to  the  dis¬ 
placements  u  and  Au  and  define  the  matrices  K0,  Ki  and  K2  as  follows 


(49) 


SifUrt  =  Y.  f  {n<1|t(u)4‘)('5u)  +  M(1)t(u)k(«u) 

m  JAm 

+Qt(u)7(Ju)}cL4  (50) 

8rjTKi(rj',  t]")t)  =  y)  f  {2  ^N*^r(u")  +  2N^1^T(u/,  u")  +  N^T(u")j  e^n^(u,^u) 

m  JAm 

+2N(1)r(£u)411)(u",  u)  +  2N(1)re(n)(u",  <5u)}  dA  (51) 

V’KafaW')*?  =  E/  {2N(n)r(<5u,u')e(n)(u",u) 

771  J  Am 

+2N(u)t((Ju,  u")eiu)(u',  <Su')}  dA  (52) 


in  which  u'  and  u"  are  allowed  to  assume  0,  u,  and  Aut-,  similarly  rf  and  Tf"  may  assume 
0,»7  i  and  A?^,  and  the  symmetric  bilinear  forms  e£n)  and  N^11)  are  defined  by 


411}(u',u")  =  ^[<MU')<MU")>  <£2(u')<Mu")>  &(u')&(u")  -t-^i(u,,)^>2(u,)]T 

N(11)  =  Aein). 


(53) 


39 


Using  the  above  definitions  in  equations  Eq.  (33)  and  Eq.  (34),  retaining  all  the  non-linear 
terms,  we  obtain  the  equilibrium  equations 


K0  +  Ki(0, 77^)  +  K  2(i?,-,»7i)  +  +  K2(t/,-,  Aifc) 


+^K2(Aj7i,A77t) 


At/,  =  AA,F 


(54) 


where  F  is  the  reference  load  vector.  The  iterative  form  of  Eq.  (54)  for  the  r-th  iteration 
of  the  i-th  step  follows: 

K;AAt7;  =  AAA^F  -  AR^  (55) 

in  which 


AAtf  =  Atf+1  -  A fft 

aaa;  =  aa:+1  -  aa; 


(56) 


and 


Kf  =  Koi  +  K1(T7i,AT/n  +  2K2(T/,,AT7n  +  K2(AT7^,AT7n 


ARr 


K 


Oi 


Koi  +  A  7,')  +  K2(,,-,  AijJ)  +  jK2(A<,  At/?) 


At/?  -  AA?F 


K0  -I-  Ki(0,  rji)  +  K2(t /,-,  T7,). 


(57) 


At  this  point,  appropriate  shape  functions  for  the  trial  and  test  functions  are  assumed 
in  order  to  discretize  equation  Eq.  (55).  In  the  quasi-conforming  shell  element  the  shape 
functions  of  Tang  [5]  are  used  for  the  discretization  of  the  linear  strain  e^,  the  total 
curvature  change  k\  the  rotations  <f>  of  the  mid-surface  and  displacements.  The  modified 
Newton-Raphson  method  and  arc  length  approach  (Riks  [4];  Zhang  and  Atluri  [6])  are 
used  to  solve  Eq.  (55).  At  bifurcation  points  on  paths,  the  asymptotic  post-buckling 
solutions  are  taken  as  the  initial  values  of  iteration,  so  that  any  assumed  imperfections 
are  rendered  unnecessary. 


5  NONLINEAR  SOLUTION  SCHEMES 

Automated  post-buckling  solution  involves:  detection  of  possible  instability  in  solution 
and  elimination  of  possible  path-retracing;  classification  of  the  detected  instabilities;  and 
computation  of  the  post-through  buckling  solution(s).  BUCKDEL  detects  solution  insta¬ 
bilities  by  monitoring  the  rank  of  the  tangent  stiffness  matrix.  Whenever  the  determinant 
of  the  tangent  stiffness  matrix  changes  its  sign  the  solution  senses  possible  instabilities  in 
that  range  of  load  and  changes  the  sign  of  the  next  load  increment  to  avoid  path-retracing. 
Through  a  cycle  of  iterations,  location  of  instabilities  are  identified  as  the  load  levels  for 
which  the  tangent  stiffness  becomes  singular.  The  identified  instability  points  are  then 
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classified  as  limit  points  or  bifurcation  points  using  some  simple  and  cost-effective  rules 
(Huang  and  Atluri  [3]).  If  the  instability  point  is  a  limit  point,  the  arc-length  continu¬ 
ation  is  enough  to  obtain  post-buckling  solution  path.  However,  if  the  instability  point 
happens  to  be  a  bifurcation  point,  the  strategies  described  in  detail  in  [3]  are  used  to 
trace  the  appropriate  post-buckling  solution  branch.  The  nonlinear  fundamental  state 
between  the  two  solution  points  n  —  1  and  n  in  the  neighborhood  of  the  bifurcation  point 
is  linearized  to  obtain  the  asymptotic  solution  for  obtaining  an  approximate  critical  buck¬ 
ling  load  factor.  A  linear  combination  of  the  normalized  eigen-vector  associated  with  the 
critical  buckling  load  factor  and  its  orthogonal  counterpart  is  used  to  used  to  determine 
the  initial  post-buckling  paths. 


6  POINTWISE  ENERGY  RELEASE  RATE  CAL¬ 
CULATION 


In  the  case  of  delamination,  the  growth  is  assumed  to  be  along  the  interlaminar  zone 
parallel  to  the  midsurface  of  the  shell  (i.e.  the  crack  cannot  shear  into  the  neighboring 
laminae).  Since  the  displacement  field  is  made  explicitly  continuous  at  the  delamination 
front,  the  delaminate  shell  cannot  slide  or  rotate  relative  to  the  base  shell.  Hence,  the 
J-integral  representing  only  self-similar  crack  growth  is  meaningful  in  the  present  case. 

The  pointwise  energy  release  rate  for  3-dimensional  crack  growth  (self-similar)  -  £/(!")) 
-  is  defined  as, 


£(r)Ar 


Limit 
e  — >•  0 


dUa 


aapni3  dx. 


dA+  f  (Ta2~^-dA  —  f  crQ2^zr~dA  (58) 

«Mi  OX  i  JA2  OX  x 


dUa 


where,  a,  fi  =  1,2, 3;  Ae  is  the  area  of  the  tube  of  radius  e  enclosing  the  crack  front;  A\ 
and  A2  are  the  areas  covering  the  ends  of  the  tube;  and  n,  a,  U  and  n  are  defined  in  the 
crack  tip  coordinate  system  x  (Figure  4). 

Consider  a  rectangular  tube  enclosing  the  delamination  front  and  passing  through  the 
nearest  stress  recovery  points  (S^)  of  the  adjoining  elements  (Figure  4).  Note  that,  the 
integrals  over  the  areas  A\  and  A2  nearly  cancel  each  other  in  a  constant  strain/stress  el¬ 
ement  model,  since  the  quantities  like  &  and  dUa/dxi  do  not  vary  in  the  neighborhood 
of  a  point  in  an  element  domain.  Then,  equation  (58)  for  the  self-similar  delamination 
growth  becomes, 

S(r)Ar  =  JArEjA  (wnx  -  aQ0hp^j  dAdT  (59) 


where  A j  (j  =  1  —  9)  axe  the  segments  forming  the  surface  area  Ae  of  the  rectangular 
tube  (Figure  4).  Since  the  Reissner-Mindlin  assumptions  (<7i3,  <j23  are  constant  over  shell 
thickness;  and  <733  =  0)  are  used  for  the  element  formulations  as  well  as  to  achieve 
displacement  continuity  at  the  delamination  edge;  and  since  n  =  {0  0  ±  1}  on  the 

segments  A  j,  j  =  4  —  9  (see  Figure  4),  the  integral  in  equation  (59)  vanishes  over 
the  segments  A  j,  j  =4  —  9.  Now,  since  n  =  {+1  0  0}  on  the  segment  Ai  and 
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n  =  {— 1  0  0}  on  the  segments  A2  and  A3  ;  equation  (59)  becomes, 

dUQ 


1  q-  J  dx$ 


dr 


£(r)Ar  =  /  [/  -/  -/  (w 

JAT  |_</  Ai  j  A2  J  A3  \ 

Now,  carrying  out  the  integration  through  the  thickness  for  each  shell,  we  get, 

S(r)Ar  =  J  rg  [w  -  (ivlQ«a>1  +  Mlaea<1  +  q13u3,i)]  dr 


(60) 


(61) 


Therefore,  as  Ar  — »  0,  the  pointwise  energy  release  rate  at  any  point  on  the  delamina¬ 
tion  front  Qg ,  computed  from  the  Gauss-point  variables,  in  a  finite  element  model  using 
constant  strain  elements,  is  given  by, 

£,(D  =Tg[W-  (NlauaA  +  MlaeaA  +  Qi3U3,i)]  (62) 


where,  JF^*)  =  (*)p(i)  -  (*)5(2)  -  (*)9{z)  and  (*)5(o  corresponds  to  the  quantities  (*) 
evaluated  at  specified  points  on  the  annular  surface.  For  example,  in  a  finite  element 
analysis,  these  specified  points  would  be  preferably  the  optimal  stress  recovery  points  - 
normally  the  Gauss  points  corresponding  to  reduced  integration  -  in  the  adjoining  element 
of  the  ith  shell  nearest  to  the  delamination  front  I\ 

Note  that  when  energy  release  rate  calculations  are  performed  with  BUCKDEL,  the 
finite  element  mesh  adjoining  the  delamination  edge  should  be  symmetric  about  the 
delamination  edge  such  that  (a)  the  stress  recovery  points  for  all  the  elements  are  at  the 
same  radial  distance  from  the  delamination  edge;  and  (b)  each  set  of  elements  from  the 
three  shells  should  have  stress  recovery  points  located  on  a  single  plane  normal  to  the 
delamination  front. 
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